/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Declarations of ray classes.

// $Id$

#ifndef __ray__
#define __ray__

#include "mcrx-types.h"

#include "blitz/array.h"
#include "mono_poly_abstract.h"


namespace mcrx {
  class ray_base;
  template <typename> class ray;
  template <typename T_numtype>
  class ray<blitz::Array<T_numtype, 1> >;
}

/** A ray_base is a base class for rays, containing basic
    functionality that has to do with positions, directions, etc, that
    does not depend on what the ray carries. */
class mcrx::ray_base {
private:
  vec3d dir_, idir_; ///< Direction inverse, precomputed to save time.
  vec3d pos_;
  T_densities n_; ///< Array of column densities of species traversed by ray.
  T_float l_; ///< Path length traversed.
  T_float d_; ///< Doppler shift factor (used when doing kinematic calculations).
  int nscat_; ///< Number of scatterings.

  void assert_stuff() {
    assert (dot (dir_, dir_) > .9999); assert (dot (dir_, dir_)< 1.0001);
  };

public:
  ray_base() {};
  ray_base (const vec3d& pos, const vec3d& dir):
    dir_ (dir), idir_ (1./dir), pos_ (pos), l_ (0), d_(1.0), nscat_ (0) {
    assert_stuff();};

  ray_base (const vec3d& pos, const vec3d& dir, const T_densities& nn):
    dir_ (dir), idir_ (1./dir), pos_ (pos), 
    n_(make_thread_local_copy(nn)), l_ (0), d_(1.0),
    nscat_ (0) {
    assert_stuff();};
  ray_base(const ray_base& rhs) : 
    dir_(rhs.dir_), idir_(rhs.idir_), pos_(rhs.pos_), 
    n_(make_thread_local_copy(rhs.n_)), l_(rhs.l_), d_(rhs.d_), nscat_(rhs.nscat_)
  {};

  /** Assignment operator. Has to make sure the column density array
      is the right size. */
  ray_base& operator= (const ray_base& rhs) {
    if (this == &rhs) return *this;
    dir_ = rhs.dir_; idir_ = rhs.idir_; pos_ = rhs.pos_; l_ = rhs.l_; 
    d_ = rhs.d_; nscat_ = rhs.nscat_;
    assign(n_, rhs.n_, false);
    return *this;
  };

  const vec3d& position () const {return pos_;};
  const vec3d& direction () const {return dir_;};
  const vec3d& inverse_direction () const {return idir_;};
  const T_densities& traversed_column() const {return n_;};
  T_float length () const {return l_;};
  T_float doppler () const {return d_;};
  /// Set the ray doppler factor to d.
  void set_doppler (T_float d) {d_=d;};
  /// Add a doppler shift of d to the ray doppler factor.
  void add_doppler (T_float d) {d_*=d;};

  int scatterings () const {return nscat_;};
  int& scatterings () {return nscat_;};

  // Note that setting the ray position resets traversed length and column.
  void set_position (const vec3d& pos) {pos_ = pos; l_ = 0; n_ = 0;};
  void set_direction (const vec3d& dir) {
    dir_ = dir; idir_ = 1./dir;
    assert_stuff();
  };
  void set_traversed_column (const T_densities& n) {
    assign (n_, n, false);};

  /** Scatters the ray in the direction defined by cos_theta and phi
      (in the ray frame). */
  void scatter_in_direction (T_float cos_theta, T_float phi);

  /// Propagates the ray a specified length, ignoring column density
  void propagate (T_float dl) {
    pos_ += dl*direction (); l_ += dl;};
  /// Propagates the ray a specified length and optical depth.
  void propagate (T_float dl, const T_densities& dn) {
    pos_ += dl*direction (); l_ += dl; 
    if(n_.size()==1) n_(0) += dn(0); else n_ += dn;};
  template<typename T>
  void propagate (T_float dl, blitz::_bz_ArrayExpr<T> expr) {
    pos_ += dl*direction (); l_ += dl; n_ += expr;};

  /** Checks if the ray intersects an axis-aligned cube between
      coordinates min and max and returns the length to
      intersection. NaN is returned if the ray misses the cube. If the
      ray is on the entry face and allow_zero is returned, 0 will be
      returned, otherwise the exit face distance will be returned. */
  T_float intersect_cube(const vec3d& min, const vec3d& max, bool assume_inside) const;
};


/** The basic ray template where the intensity is a scalar. */
template <typename lambda_type>
class mcrx::ray: public ray_base {
public:
  typedef lambda_type T_lambda;
private: 
  T_lambda i_; ///< The intensity of the ray
public:
  ray() {};
  ray (const vec3d& p, const vec3d& d, const T_densities& nn, 
       const T_lambda& ii):
    ray_base (p, d, nn), i_ (ii) {};
  ray (const vec3d& p, const vec3d& d, const T_lambda& ii):
    ray_base (p, d), i_ (ii) {};
  const T_lambda& intensity () const {return i_;};
  void set_intensity (T_lambda i) {assign(i_, i);};
  void bias(const T_lambda& factor) {
    assert(condition_all(factor>=0));
    assert(condition_all(factor<blitz::huge(T_float())));
    ;i_ *= factor;};
};


/** Specialization of the ray template for the polychromatic case when
    the intensity is a blitz array. */
template <typename T_numtype>
class mcrx::ray<blitz::Array<T_numtype, 1> >: public mcrx::ray_base {
public:
  typedef typename blitz::Array<T_numtype, 1> T_lambda;
private: 
  T_lambda i_; ///< The intensity of the ray
public:
  // the assignment operator is fully functional so we can have a
  // default constructor
  ray () {};
  ray (const vec3d& p, const vec3d& d, const T_densities& nn, 
       const T_lambda& ii):
    ray_base (p, d, nn), i_ (make_thread_local_copy(ii)) {};
  ray (const vec3d& p, const vec3d& d, const T_lambda& ii):
    ray_base (p, d), i_(make_thread_local_copy(ii)) {};
  ray (const ray& rhs) : ray_base(rhs), 
			 i_(make_thread_local_copy(rhs.intensity())) {};

  /** Assignment operator. Has to make sure the intensity array is the
      right size. */
  ray& operator= (const ray& rhs) {
    if (this == &rhs) return *this;
    ray_base::operator= (rhs);
    assign(i_, rhs.i_, false);
    return *this;
  };

  const T_lambda& intensity () const {return i_;};
  void set_intensity (T_float i) {i_ = i;};
  void set_intensity (const T_lambda& i) {assign(i_, i);};
  template<typename T>
  void set_intensity (blitz::_bz_ArrayExpr<T> i) {i_= i;};

  void bias(T_float b) {
    assert(b>=0);
    assert(b<blitz::huge(T_float()));
    i_ *= b;};
  void bias(const T_lambda& b) {
    assert(condition_all(b >= 0));
    assert(condition_all(b < blitz::huge(T_float())));
    i_ *= b;};
  template<typename T>
  void bias(blitz::_bz_ArrayExpr<T> b) {
    assert(condition_all(b >= 0));
    assert(condition_all(b < blitz::huge(T_float())));
    i_ *= b;};
};

#endif
